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Abstract. Motion of stellar-mass satellites is studied around a massive compact body which is surrounded by a 
gaseous slab of a stationary accretion disc. The satellites suffer an orbital decay due to hydrodynamical interaction 
with the disc medium (transitions across the disc, gap opening in the disc, density waves) and gravitational 
radiation. Arbitrary orbital eccentricities and inclinations are considered, and it is observed how the competing 
effects depend on the parameters of the model, namely, the mass and compactness of the orbiters, the osculating 
elements of their trajectories, and surface density of the disc. These effects have a visible impact on the satellites 
long-term motion, and they can produce observational consequences with respect to galactic central clusters. It 
is shown that the satellite-disc collisions do not impose serious restrictions on the results of gravitational wave 
experiments if the disc medium is diluted and the orbiter is compact but they are important in the case of 
environments with relatively high density. We thus concentrate on application to accretion flows in which the 
density is not negligible. We discuss the expected quasi-stationary structure of the cluster that is established on 
sub-parsec scales within the sphere of gravitational influence of the central object. Relevant to this region, we give 
the power-law slopes defining the radial profile of modified clusters and we show that their values are determined 
by satellite interaction with the accretion flow rather than their initial distribution. 
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1. Introduction 

There is ample evidence that central regions of galaxies 
consist of three main constituent parts: a supermassive 
dark object (presumably a black hole), a dense cluster of 
stars, and an accretion flow with disc-type geometry (for 
a review, see Rees 1998; Kato, Fukue & Mineshige 1998). 
Related to these three ingredients are typical length- 
scales: (i) r g = GM/c 2 w 1.5 x 10 13 M 8 cm (gravita- 
tional radius of the central black hole expressed in terms 
of its mass M 8 =M/1O 8 M < 1). (ii) r c = GM/a 2 w 
0.5Af§cr^QQQ pc, the size of an inner star-cluster (where 
motion of its members is governed mainly by the black- 
hole gravitational field) written in terms of stellar velocity 
dispersion criooo = cr/(1000km/s). A typical number of 
stars in this region can be of the order N « 10 5 . Finally, 
(hi) rd « 10 4 T g (the size of the disc). For the purpose of 
this paper these stars play a role of stellar-mass satellites 
orbiting the center. The disc itself will be described in the 
thin-disc approximation with the total mass Md<C-M. 

In this paper we assume that the three mentioned sub- 
systems coexist and they form an integral structure for a 
sufficiently long period of time. This is a reasonable as- 
sumption for a range of model parameters (discussed be- 
low), and the resulting geometry captures interesting ef- 



fects of realistic systems, although it is simplified in several 
respects. The three subsystems are of very different na- 
ture, so that their symbiosis and mutual interactions pose 
a complex, ecological problem (e.g., Shlosman & Noguchi 
1993; van der Marel & van den Bosch 1998; Zwart, Hut 
& Verbunt 1997). We confine ourselves to a simple analy- 
sis in which the gravitational field of the system is dom- 
inated by the central compact mass while long-term dy- 
namics of the satellites is influenced by dissipation in the 
disc medium and by gravitational radiation. The present 
paper is motivated mainly by indications of substantial 
concentrations of both stars and interstellar medium in 
the galactic central parsec. For example, a toroidal sys- 
tem of stars and gas has been proposed as a plausible 
explanation of brightness peaks observed in the nucleus of 
M31 (Tremaine 1995; Bckki 2000). However, as we want 
to discuss also the regime in which inspiraling satellites 
are influenced by gravitational radiation in the field of 
the massive centre, we need to consider sub-parsec scales 
which are not yet resolved observationally. This appears 
relevant for the continuation (towards small length-scales) 
of the velocity dispersion versus mass relation (Ferrarese 
& Merritt 2000), which is so important in the discussion 
of the masses of supermassive black holes. Finally, core- 
collapsed globular clusters could be another class of rele- 
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vant objects on scales of masses w 10 4 Mq, i.e. much less 
than galactic cores. 

We proceed along the same course as several recent 
works in which further references can be found: Rauch 
(1999) studied stellar dynamics near a massive black hole, 
including the effects of general relativity on stellar trajec- 
tories; Subr & Karas (1999) examined the long-term evo- 
lution of orbital parameters of a satellite colliding with a 
thin Kcplcrian disc; Narayan (2000) estimated the relative 
importance of hydrodynamic drag versus gravitational- 
radiation decay of the satellite orbit (this author was 
interested primarily in the case of a body embedded in 
an advection-dominatcd flow very near the central hole); 
Collin & Zahn (1999), on the other hand, explored the fate 
of satellites moving on the periphery of the disc and at in- 
termediate distances (rs 1 pc for Mg w 1) where mutual 
interaction of the disc matter with the satellites triggers 
the starburst phenomenon and provides a mechanism of 
metal enrichment; finally, Takcuchi, Miyama & Lin (1996) 
and Ward (1997) discussed the problem of satellite migra- 
tion via gap formation and density waves in the disc, and 
they gave a detailed overview of related topics which have 
been widely discussed also by others. 

The orbital motion of satellites forming the central 
cluster is dominated by gravity of the central body; how- 
ever, the long-term evolution of their trajectories is af- 
fected also by many other effects such as gravity of the 
background galaxy (e.g., Sridhar & Touma 1999), dy- 
namical friction and tidal interactions (Colpi, Mayer & 
Governato 1999), self-gravity of the disc matter (Shlosman 
& Begelman 1989; Vokrouhlicky & Karas 1998), repetitive 
collisions with an accretion disc (Syer, Clarke & Rees 1991; 
Rauch 1995), and gravitational radiation losses (Peters 
& Mathews 1963). The interaction between the members 
of the cluster and the disc are of interest also because 
they may gradually change the composition of the disc 
(Artymowicz, Lin & Wampler 1993), speed up the evolu- 
tion of the satellite itself (Collin & Zahn 1999), and they 
can enrich the environment outside the disc plane (Zurck, 
Sicmiginowska & Colgate 1994; Armitage, Zurek & Davics 
1996). 

Using simple estimates we compare the energy and an- 
gular momentum losses via gravitational radiation against 
the hydrodynamical drag acting on the satellites. We as- 
sume that their orbits have, initially, arbitrary eccentrici- 
ties and nonzero inclinations with respect to the disc plane 
- the situation which is complementary to the case of 
planetary formation and their migration inside the disc. 
For the subsequent (low-inclination) stages of the orbital 
evolution we employ a simple prescription in which grav- 
itational radiation losses still tend to bring the orbiter 
towards the centre while hydrodynamical effects are ap- 
proximated in terms of motion through gas. 

We build our discussion on a model description (Subr 
& Karas 1999) in which the satellite suffers from dissipa- 
tion of the orbital energy during repetitionary star-disc 
interactions but the disc itself does not change with time. 
We show evolutionary tracks followed by the satellite in 



the parameter space of its orbital osculating elements. The 
results are found quite sensitive to actual values of the 
disc density, the mass of the satellite, and to several other 
parameters of the model. We discuss the dynamical evolu- 
tion of a cluster of satellites which gradually departs from 
its initially spherical shape and evolves into a flattened 
system. The initial stage of the evolution (when most of 
the orbits arc still inclined with respect to the disc) is 
similar to the situation discussed by Pincault & Landry 
(1994) and Rauch (1995), but for the later stages we must 
adapt the prescription for star-disc interactions. See also 
Ivanov, Papaloizou & Polnarev (1999) for a complemen- 
tary discussion of the gradual evolution of a circumbinary 
disc structure. 

2. Satellite motion influenced by the disc 

The satellite body is assumed to cross the disc at super- 
sonic velocity with a high Mach number, M. « 10 2 -10 3 . 
We recall (Syer et al. 1991; Vokrouhlicky & Karas 1993; 
Zurck et al. 1994) that passages last a small fraction of 
the orbital period at the corresponding radius where they 
occur, and they can be considered as instantaneous, repet- 
itive events at which the passing satellite expels from the 
disc some amount of material that lies along its trajectory. 
In terms of inclination angle i and the disc thickness h, 
the typical ratio of the two periods is 6~h/(rsmi); this 
quantity is assumed to be less than unity. Hence, for a 
geometrically thin disc (e = h/r <C 1), its influence upon 
inclined stellar orbits can be treated as tiny kicks (im- 
pulsive changes of their energy and momentum) at the 
points of intersection with the plane of the disc. This can 
be translated using the relation for the speed of sound, 
c s = eilr, to the claim that the motion across the disc is 
indeed supersonic, which is also required for consistency. 
Smallncss of S means that the satellite remains outside the 
disc for most of its revolution around the centre; naturally 
this condition cannot be ensured at final stages when i is 
very small. 

2.1. Collisions with the disc 

Repetitive collisions lead to the gradual change of the or- 
biter's velocity v — ► v' , which can be expressed by mo- 
mentum conservation (Subr & Karas 1999): 

(A + 1) v' = v r e r + v$e# + (v v + Av K ) e^ 



(1) 



in spherical coordinates (■& — ir/2 is the disc plane). Here, 
A(r) = EdfreiS" 1 ^ , w rc i is the relative speed between 
the orbiter and the disc matter, £* is the column density 
characterizing the compactness of the orbiter and defined 
by £* = M*/ (ttR-I) (quantities denoted by an asterisk 
refer to the orbiter, M r ^_M), and Ed(r) is the disc surface 
density. Rotation of the disc is assumed Keplerian, v^, = 

Let us consider a satellite on an orbit with semi-major 
axis a, eccentricity e, inclination i, and longitude u) of the 
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ascending node. Eq. (uh implies a set of equations which 
can be solved numerically in terms of the orbiter's oscu- 
lating elements, while analytical solutions are possible in 
special cases (Subr & Karas 1999). As a useful example, 
we assume a power-law surface-density distribution in the 
form Ed = K(r/r s ) £q (K — const) and we adopt a 
perpendicular orientation of the orbit, cos a; = 0. We find 



dx+ (x\ + C 2 xJ) , y=l + C 2 x. 



-3 



(2) 



where x± = l±x, x = cos A y = 1 — e 2 . Strictly speaking, 
this formula concerns only the case of orbits intersecting 
the disc at two points with identical radial distances from 
the center but it can be used also as an approximation 
for orbits with arbitrary orientation. We remark that the 
relative accuracy Aa/a of the determination of the semi- 
major axis is better than 15% with reasonable density 
profiles (s of the order of unity) ; numerical computations 
are not limited by assumptions about u> imposed in eq. (0). 
Furthermore, one can write 



KHr. 







3 — y — 2x 



1 



(3) 



Integration constants C\^ are to be determined from 
initial values of a = ao, x = xq, and y = yo. Then, the 
temporal history is obtained by integrating over the or- 
bital period, dt = 2na 3 / 2 />/GM, in the form 



2tt 



a(x) 3/2 [a(x)y(x)/r g ] s dx 
KVGM S© J Xo ^{3-y(x)-2x)(l-x 2 ) 



(4) 



Here, a factor missing in eq. (19) of Subr & Karas (1999) 
is corrected (no other equations and graphs were affected 
by that omission) . The orbital decay manifests itself in the 
gradual decrease of a, e and i, for which surface density 
of the disc is the main factor. The time derivative of the 
semi-major axis is 



flcoi = By 



where 



-in 



^0 



'3 — y — 2x .„ . 

,, 2 . (2-x-y), 
2/(1 -* ) 



B=~BcM§ 1 fi q E 2 



0.1 



</:f 



(5) 



(6) 



and fi E = M/(0.1Me) is the accretion rate in units of 
Eddington accretion rate (with a 10% efficiency factor in- 
troduced). The factor B and power-law indices qi..A are 
determined by details of the particular model adopted to 
quantify the disc properties. Table ffl gives the values rele- 
vant for different regions of the standard disc as well as for 
the gravitationally unstable outer region (Collin & Hure 
1999). Notice that the algebraic functional form of radial 
dependencies remains identical in all these cases (a power- 
law), and we can use it with convenience also later for 
different prescriptions of the satellite-disc encounters. 



Table 1. Parameters in eqs. (||), (0) and (|9|) describing 
the orbital decay in the case of different regimes and for 
different disc models. 
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Notation used in table. 

Disc models: (i) ... Standard disc with p = p la ,d, s — 3/2; 
(ii) . . . Standard disc p = p gas , S = —3/5, electron 
scattering opacity; (iii) . . . The same as (ii) but with 
s = —3/4 and free-free opacity; (iv) . . . Marginally unsta- 
ble self-gravitating disc (solar metallicity, optically thick), 
s — —15/7; (v) . . . The same as (iv) but for zero metallic- 
ity, optically thin medium, s = —5/3; cf. Sec. H for further 
details. 

Regimes of orbital decay: Orbital decay dominated by star- 
disc collisions (col), by gap formation in the disc (gap), 
and by density waves (dw), respectively. 



2.2. An orbiter embedded in the disc 

As drag is exerted on the satellite body, its orbit becomes 
circular and declined in the disc plane. Within this frame- 
work, orbital eccentricity and inclination are expected to 
reach zero values in the disc, so that quasi-circular trajec- 
tories are relevant near the center. 

The orbit evolution is thus reduced to the situation 
which was addressed by several people (e.g., King & 
Done 1993; Takeuchi et al. 1996; Ward 1997; Ostriker 
1999) in connection with formation and subsequent mi- 
gration of bodies inside the disc. Two basic modes 
can be distinguished according to the disc properties 
and the orbiter mass. First, a gap is cleared in the 
disc if the satellite's Roche radius exceeds the disc 
thickness, r L w (M*/M) 1 / 3 r^/i, and simultaneously 
M„>Af gap wV40ae 5/2 M (Lin & Papaloizou 1986). Motion 
of the satellite is then coupled with the disc inflow, so that 



B. 



(7) 



On the other hand, if the satellite is unable to create the 
gap, the gas drag is imposed on it through quasi-spherical 
accretion. The resulting radial drift is weaker by a factor 
e <C 1 than the drift caused by density-wave excitation 
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Fig. 1. Contours of / are plotted as a function of incli- 
nation and eccentricity for the case (ii). Function f(i,e) 
characterizes the relative importance of energy losses in 
eq. (flq). Values of / are indicated with contour lines. A 
saddle-type point develops between two neighbouring con- 
tours, / = 0.150 and / = 0.152, which are plotted with 
dashed lines. Its exact location depends on the disc model 
but the overall picture remains very similar in other cases, 
too. 

(Ward 1986; Artymowicz 1994). Hence, such a satellite 
migrates inward mainly due to the latter effect on the 
time-scale 



tdv 



{CM*n) 



1 ,2 



M 



f2 



7rr 2 E r 



(8) 



where C is a dimensionless constant of the order of unity. 
Using eq. (H) and da/dtmr/tdw we obtain 



M* 



«dv 



M, 



B. 







(9) 



Substantial differences in the satellite migration are thus 
introduced in the model already within this very simpli- 
fied picture where the process of satellite sinking is driven 
by the gas medium. For the region of the gas pressure 
dominated standard disc and for a self-gravitating zero- 
metalicity model (v), the chance of opening the gap in- 
creases with decreasing radius. The situation is opposite 
in the case (i) (radiation pressure dominated) and for the 
solar-metallicity disc (iv). 

2.3. The orbital decay due to gravitational radiation 

The orbiting companion emits continuous gravitational 
radiation whose waveforms are of particular relevance 
for gravitational wave searches from compact binaries 
in the Milky Way. Possible approaches to their observa- 
tional exploration have been discussed by several people 
(e.g., Nakamura, Oohara & Kojima 1987; Dhurandhar & 
Vecchio 2001; Hughes 2001). The average rate of energy 
loss which the orbiter experiences via gravitational radia- 



tion over one revolution can be written in terms of orbital 
parameters (Peters & Mathews 1963), 

. 3 2 G « M ne - 1 + ge2+ g e . ) . (10) 



E — 



5 C 5 a 5 (! _ e 2)7/2 

Corresponding to eq. (fill) are the change of semi-major 
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The above formulae (|TC])-(|l^) assume that the satellite 
star follows an eccentric orbit in Schwarzschild geome- 
try of the central massive body. Gravitational radiation 
losses compete with those caused by star-disc encoun- 
ters. We are thus interested in the relative importance 
of these mechanisms, which is characterized by the ratio 
ftcol/gw = [da/di] col / [da/dt] gw : 



K, 



1/gv 



col/gw 



hK M 
32^ Ml 



£fr 







-I s+5/2 



f(x,e), 



(13) 



where, for whichever of the models described by the power- 
law density profile, 



f(x,e) 



{l + e 2 -x)(l-e 2 Y +5/2 J2 + e 2 -2x 
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(14) 



For the standard thin disc model (i) we obtain 



K {i \, « 1.8 x lO^MguZ 1 

col/gw °"e 
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-l 
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— f(x,e), 
while for a gas pressure dominated disc (ii) 
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For the model (iii) we find 
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/(x,e). (16) 



f{x,e) 



(17) 



We show the functional form ( |14|) in Fig. [j] where the fac- 
tor f(x,e) determining the dependence of 7\l co i/ gw on incli- 
nation and eccentricity is plotted for the case (|l^) . In other 
words, /(x, e) represents that part of the drag ratio that 
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is independent of the disc medium and the satellite physi- 
cal properties; only the two mentioned orbital parameters 
play a role here. Typically, for a solar-type star it is only 
on eccentric orbits that / becomes small enough to bring 
7^-coi/gw below unity. The required eccentricity is rather 
high, and such a satellite would be trapped or disrupted 



directly by the central hole. Otherwise, 1Z 



col/gw 



> 1 



for a^r g and for standard values of the disc parameters 
(a Si 1, ti E ~ 1). This means that direct hydrodynam- 
ical interaction with the disc plays a dominant role in 
the orbital evolution of satellites crossing the disc, unless 
the medium is extremely rarefied (e.g., an ADAF; Narayan 
2000). Notice that the point /(l, 0) = (i.e. a fully circu- 
larized orbit inclined into the disc plane) is the exception, 
where the adopted approximation of instantaneous colli- 
sions breaks down. 

Analogous to 7?. co i/ gw , one could explore the relative 
ratio of the hydrodynamical versus gravitational radiation 
losses in other regimes of the satellite-disc interaction. In 
this way, the relevant formulae (those which apply in the 
course of orbiter evolution) are 
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for the inner, radiation pressure dominated disc (i), and 



^dw/gw = 1-0 x 10~ 
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, (20) 
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for the middle region (ii). Here, the eccentricity-dependent 
factor was omitted upon the finding that orbits are almost 
circular when inclined in the plane of the disc. It is evident 
from cqs. (p3[) (]21[) that gravitational radiation can have 
a visible impact only at small a, especially when eq. ( |l8| ) 
applies (cp. also Fig. Sand related discussion below). 




3x10" 



Fig. 2. Evolutionary tracks of the satellite are plotted 
in the plane semimajor axis a (units of r g ) versus time 
(years). Two very different cases are shown, (a) Left: solar- 
type satellite (£* = Eq). The three curves correspond 
to different initial eccentricities (from top to bottom), 
e = 0, 0.4, and 0.8. Each curve is divided into three seg- 
ments. Solid line indicates the period of motion outside the 
disc (when collisions occur with large orbital inclination) ; 
dashed line corresponds to density-wave driven motion in 
the disc; dotted line is the late stage with gap formation, 
(b) Right: a compact satellite (£» = 10 7 Eq). Only the 
first stage is resolved here. In all cases, ao — 100, xq = 0. 
Gravitational radiation contributes to orbital changes in- 
dependent of inclination but its impact is clearly impor- 
tant near the centre, where a decreases rapidly. 



neglected. The satellite sinks in the disc where the impul- 
sive approximation (H) loses its validity and it is substi- 
tuted by motion in the disc plane. Time-scales are gener- 
ally longer in the latter, compact satellite case, although 
gravitational wave emission speeds up the evolution at 
very late stages (a ^ 50). For definiteness, we adopt the 
disc model (iii) and the condition of tan« tr = h/r when 
the transition in the disc occurs. We checked by modifying 
the value of tani t r by a factor of 10 that the qualitative 
picture of orbital evolution does not depend on its exact 
choice and also that numerical results remain similar. 



2.4. Evolutionary tracks of the satellite 

Now we explore the evolutionary tracks of the satellite in 
the parameter space of osculating elements. We start by 
considering two effects: gravitational- wave losses in the ap- 
proximation of eqs. (JK^)-(|l2|), and hydrodynamical drag 
acting on the satellite according to eq. (|l|) twice per revo- 
lution. Dissipation operates with an efficiency depending 
on the type of satellite, and it provides a mechanism for 
the separation of different types of bodies in the phase 
space of a cluster. 

Typical results of orbit integrations are presented in 
Figure 0. Notice the big difference in time-scales relevant 
for non-compact stars (left panel) when compared with 
compact ones (right panel). In the former case, hydrody- 
namical drag is more pronounced. It gradually changes the 
orbital plane, while gravitational radiation can be safely 



2.5. Limitations of the adopted approximations 

For simplification, we ignored the effects of gradual change 
of the mass of both the satellite and the central body. 
Furthermore, we assumed that the interaction has no ef- 
fect on the structure of the satellites. This is a plausible 
assumption for stars with column densities much larger 
than that of the disc, while it is inadequate for giants 
which must quickly lose their atmospheres (Armitage et 
al. 1996). Also, we did not consider various effects act- 
ing on the disc structure (e.g. torques imposed on it by 
the dense cluster of stars; Ostriker 1983). Although all 
these effects will be important for a complete unified treat- 
ment of accreting black holes in active galactic nuclei, they 
can be neglected without losing the main physical effects 
influencing the satellite motion in the present simplified 
scheme. 
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We note that the significance of direct orbiter's col- 
lisions with the disc material is controlled by a charac- 
teristic time-scale r oc £*/£<} (r), which we expressed in 
terms of the orbital decay a. Let us recall that dimcn- 
sionlcss parameter characterizing compactness of the or- 
biter can be introduced in different ways. While £*/£d 
stands directly in the description of star-disc collisions, 



the usual factor e — GM*iL 



determines the impor- 



tance of general relativity effects near the surface of a 
compact body: e w 1CT 6 (S*/S ) (R*/Rq). We treat 
motion in the Newtonian regime and we only take the 
possibility of the satellite capture into account by remov- 
ing the orbiter from a sample if its trajectory plunges too 
close to the central mass, below a marginally stable orbit. 
Another dimcnsionless quantity has also been designated 
as the compactness parameter when describing accretion 
onto compact objects, e = Lcr T r _1 ?ri ( r 1 c _3 . It considers 
the effect of radiation luminosity L acting through a cross- 
section <7 T in the medium, however, we can safely ignore 
radiation pressure on macroscopic satellites hereafter. 

The gradual and monotonic decrease of eccentricity 
is overlaid with short-term oscillations if the disc mass 
is non- negligible (Vokrouhlicky & Karas 1998). Also the 
satellites' inclination converges to a somewhat different 
distribution (instead of a strictly flattened disc-type sys- 
tem) when two-body gravitational relaxation is taken into 
account (Subr 2001; Vilkoviskiy 2001 - preprint). On the 
other hand, complementary to the scenario of the satel- 
lites grinding into the disc is the picture of enhanced star 
formation in the disc plane, in which case the stars are 
born with zero inclination (Collin & Zahn 1999). But these 
effects, as well as evaporation processes operating in the 
cluster, as suggested by various Fokker-Planck simulations 
(e.g., Kim, Morris & Lee 1999), remain beyond the scope 
of the present paper. 

We could see that different mechanisms (of which we 
considered particular examples) affect the orbital motion 
rather selectively, depending on the orbiter's size and 
mass. One thus expects separation of different objects in 
the cluster phase space. In order to verify this expectation 
we examine in the following paragraph a simple scheme 
which captures gradual changes in the structure of the 
cluster. Such a discussion is required: indeed, in the ab- 
sence of sufficient resolution which would enable tracking 
of individual stellar paths in nuclei of other galaxies, one 
needs to inspect the overall influence on the members of 
the cluster, namely, the change of the radial distribution 
of the satellites in terms of their number fraction and av- 
erage inclination. 

Let us note that the accretion flow is supposed to re- 
main undisturbed by the presence of the embedded clus- 
ter. This assumption gives an upper limit on the total 
number of stars inside the radius r^ (s=s 10 4 r g ), and on the 
fraction of those dragged into the disc plane in this region. 
A simple smooth disc can be destroyed, especially in the 
process of gap formation (the case of sufficiently large M* 
and small h); the models (i) and (v) are particularly sus- 
ceptible to the occurence of multiple gaps. Very roughly, if 



~ 2% of the total number of N — 10 4 satellites get aligned 
with the disc at late times (a result of our computations 
for the model (iii)), then their Roche lobes might cover an 
area of the order ss 0.01A?>dr L . This is just comparable 
with the total disc surface for M* ss IMq, however, re- 
call that there are more conditions for the gap formation 
depending on the disc model. Also, effects of two-body 
relaxation and of satellite scattering in the gravitational 
potential of the disc are expected to decrease the fraction 
of aligned bodies when these effects are taken into account. 

3. The cluster evolution 

We follow the cluster evolution starting from a clean sta- 
tionary state (no dissipative perturbation), which corre- 
sponds to a spheroidal stellar system gravitationally cou- 
pled with the central black hole (Bahcall & Wolf 1976). 
Hence, number density no(a) of the satellites is assumed 
proportional to a 1 ' 4 , while the initial distribution in ec- 
centricities and inclinations conforms with the condition 
of a randomly generated sample. Common parameters of 
the computations are the mass M = 10 8 Mq of the central 
body, and the orbiters' column density £* = £q. Now, 
the perturbation effects are switched on, depending on the 
prescribed disc type and accretion rate (typically /i E = 1 , 
but we also examined other situations and less conven- 
tional values of the above parameters; cf. Subr 2001). 

As a consequence of the dissipative action, the system 
departs from the initial configuration. We concentrate our 
main attention on the range of radii r Ss 10 3 r g , well below 
the outer edge of the disc r^ w 10 4 r g (where Ed becomes 
negligible) . We find a new quasi-stationary, modified clus- 
ter distribution developing in that region. On the outer 
boundary, fresh satellites could be inserted in the system 
from an external reservoir in order to maintain the steady 
flow towards the center, but in this work we alternatively 
consider a closed system with a large but finite number 
of its members, Nq, which are not replenished. The ac- 
tual value of iVo is not very important in the present work 
because we neglect gravitational interaction among the 
satellites themselves (No will play a role when two-body 
collisions are taken into account). 

3.1. A role of the disc model 

Several different models of the disc were adopted. We spec- 
ify them in terms of their corresponding density profiles, 
£d(r), and vertical thicknesses, h(r): 

(i) Radiation pressure dominated, geometrically thin disc 
(Shakura-Syunaev a-viscosity model; Frank et al. 1992) 
(ii) Gas pressure dominated standard disc with opacity 
due to electron scattering. 

(iii) The same as in (ii) but with opacity dominated by 
free-free scattering. 

(iv) Gravitationally unstable discs in the region of a self- 
gravitating disc beyond the Toomre radius; see Sec. (2.3) 
of Collin & Hure (1999) for analysis of the marginally un- 
stable (£ = 5) self-gravitating disc with solar metallicity 
and large optical thickness. 
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Fig. 3. Number densities nt(a) of the cluster. Left panel: 
quasi-stationary states for M* = IMq (solid line), 3Mq 
(dashed), and IOMq (dotted). Right panel: snapshots of 
n t (a) are shown for M* = IMq at f = 10 6 (solid), 10 8 
(dashed), 10 10 (dotted), and 10 12 (dash-dotted). In all 
cases, £* = £q (solar-type compactness of the satellite). 
On the vertical axis, units are arbitrary and scaled with 
respect to the initial profile no (a) ex a 1 ' 4 . Units of GM/c 2 
are used for a on the horizontal axis. 



(v) The same as in (iv) but for optically thin medium 
with zero metallicity. 

The choice of different models is dictated by their appli- 
cability for the description of accretion regimes relevant 
under different conditions, namely, the cases (iv)-(v) are 
appropriate for black-hole discs at large (<; 10 3 r g ) radii. 
One could consider other models relevant in various situa- 
tions but the adopted examples cover qualitatively differ- 
ent evolutions which are frequently encountered. In partic- 
ular, the above models comprise a range of Ed/S* depen- 
dencies (dimensionless factor determining the evolution- 
ary time-scales) and different (/i/r)(Sa/S*) = e/V whose 
radial dependence determines area in the disc where the 
gap is formed. 

First, we carried out simulations while sticking to one 
of these models in the whole range of radii from the centre 
up to the outer edge, even if it extends somewhat out of 
the zone of validity of the particular model. This helps us 
to avoid additional complexities which could stem from a 
complicated prescription for the disc properties. For ex- 
ample, one expects a different rate of radial drift across 
boundaries where the disc properties are switched, but we 
refrain from this complication for the moment. This ef- 
fect will be observed later on: to achieve a more realistic 
description, a further step will be carried out by joining 
different regions of the disc in a unified model. In that case, 
the radiation pressure supported standard disc might be 
limited to the innermost part of the system while gravita- 
tionally unstable part takes over at distances of the order 
of a few hundreds r g . 

Parameters defining radial migration in the models (i)— 
(v) were given in Table [fl. Here we only remark that the 
factor B determines the pace in absolute time units with 
which the evolution proceeds, while power-law indices qi 
influence the form of the satellites number density distri- 
bution n(o). 



3.2. Modified cluster distributions 

The main results are shown in Figure g, in which the disc 
model (iii) is adopted. This will serve as a representa- 
tive test case (more examples are deferred to Appendix). 
It is particularly interesting to observe transition regions 
where the mode of the mutual satellite-disc interaction is 
changed. At the transition, the efficiency of radial transfer 
of the satellites is changed, and this effect produces sudden 
changes (either dips or concentrations) in n t {a) profiles. If 
the average inclination of the satellites is small (the case 
of a flattened configuration) , a ring is created in the place 
where the orbital decay is stalled. 

The computations were launched from the initial dis- 
tribution no(a) ex a 1 ' 4 . In the left panel, Fig. Ha, new 
nt(a) profiles are shown at a sufficiently late time, t = 
10 12 , when a quasi-stationary state has been established.FI 
Three cases are plotted with increasing satellite masses: 
M* = IMq, 3Mq, and IOMq. After being inclined to the 
disc, low mass satellites proceed via density waves (e.g., 
1M* case with n t (a) oc a" 1 / 2 at small a). On the other 
hand, very massive satellites develop a gap already very 
far from the centre and they proceed in this mode the 
whole way down (10M* case with nt(a) oc a 1 ' 4 , which 
is incidentally parallel with the initial no distribution).^ 
Intermediate satellites show an evident dip in n t {a): its 
formation is connected with the gap in the disc which 
arises at sufficiently small a. Exact location of the dip de- 
pends on the disc model too; for 3Mq satellites it occurs 
at a < 150. 

The right panel, Fig. |^b, is complementary to the left 
one, showing temporal evolution of n t (a) for one solar 
mass satellites. One notices two areas with different rit(a) 
slopes in the graph: at small a the overall distribution is 
eventually dominated by satellites residing in the disc (cor- 
responding to negative slope), while a whole mixture of 
satellite inclinations persists farther out (a positive slope). 
Transition between the two regions occurs at a ~ 10 2 (its 
location moves gradually towards larger a as time pro- 
ceeds; see the figure). A terminal dip is seen at large a (on 
the right of the graphs); it is caused by depletion of the 
sample in the simulation. 

A realistic situation must involve further changes of 
the nt(a) slope caused by varying the accretion mode 
(see Appendix and Fig. |A.l| for details). In this respect, 
the form of quasi-stationary distributions is of particular 
interest (the persisting quasi-stationary profile of n t {a) 
is given more detailed discussion in the next section, 
Sec. ^^) . We remark that the long-term evolution of the 



1 Time is given in units GM — c — 1, total number of satel- 
lites is N — J rit (a) da. Conversion of time t to physical units 
is achieved by i p hys[yr] = 1.6 x 10 _5 A/gi, while for distances 
ttphysfcm] = 1.5 x 10 13 A'/ 8 a- 



The values of nt{a) slopes are further explained in Sec. 3.3 



We recall that these slopes refer to the motion in the disc. The 
overall nt(a) profiles can be thus different in the region where 
the whole sample is dominated by those satellites moving along 
inclined trajectories outside the disc. 
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Fig. 4. Left panel: the shape of the distribution n t (a) 
of a cluster containing compact orbiters of the masses 
M* = 1.5Mq (solid line), 3M Q (dashed) and 10M Q (dot- 
ted) at t = 10 12 . Right panel: subsequent phases of evolu- 
tion for a cluster of neutron stars with M* = 1.5Mq. The 
power-law is represented by a line in the log-log plane of 
Tit (normalized by iVo, the initial number of cluster mem- 
bers) vs. a (units of r g ). Again, a thin dotted line stands 
for the initial power-law distribution, Hq. Other curves cor- 
respond, from top to bottom, to t = 10 9 , 10 11 , 10 13 , 10 15 , 
and 10 17 . 
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Fig. 5. Left panel: a quasi-stationary n t (a) profile is plot- 
ted (by crosses) for the model (i), t = 10 13 . A power-law 
approximation, nt(a) oc a qi , is constructed in appropri- 
ate ranges of radius (solid lines). Suppression of the profile 
occurs at small a in consequence of quite efficient gravi- 
tational radiation losses near the center. Right panel: the 
slope evolves from the initial value, q — 1/4, until the 
quasi-stationary state of q — > q{ is reached. Two snapshots 
are shown for the model (iii): t = 10 10 (upper profile), and 
t = 10 12 (lower profile). See the text for corresponding val- 
ues of <Jf . 



orbits is influenced by collisions with the disc especially 
if the satellite's cross-sectional area is large, e.g. a giant 
star or a solar-size body, but quite the opposite conclusion 
can be drawn for compact objects such as neutron stars 
and stellar-mass black holes for which the cross-section 
is much smaller, diminishing their collisional interaction 
with the disc. The whole cluster is modified by this kind 
of dissipative process in a selective manner, depending on 
M*. 

In Figure [|a we present the resulting form of the distri- 
bution n t —iQi2 (a) for different masses of compact bodies. 
On the other hand, panel [|b concerns a cluster consisting 
of one solar-mass satellites, M* = IMq, shown at differ- 
ent times (several consecutive moments are examined with 
a logarithmic spacing in t). Now, there is no need to con- 
sider separately that part of the distribution of trajecto- 
ries inclined into the disc; indeed, no compact bodies reach 
the disc plane, even at late phases of the evolution. This 
implies that the orbits evolve predominantly by gravita- 
tional radiation and the resulting graphs come out almost 
identical for all the disc models under consideration. We 
are thus led to restrict ourselves to non-compact satellites 
which suffer relatively intense drag from the disc. 

3.3. A quasi-stationary state 

Time-scales of the perturbed cluster evolution depend on 
the disc model and, as shown above, a number of other 
factors characterizing the cluster members - their typical 
mass, compactness and efficiency of the dissipative inter- 
action with the disc medium. With solar-type (and less 
compact) satellites, the disc plays a dominant role in de- 
termining the form of rit(a). The starting canonical distri- 
bution no(a) oc a 1 / 4 (corresponding to p oc r~ 7 / 4 in the 



notation of Bahcall & Wolf 1976) changes to nt(a), where 
it stands almost frozen in a quasi-stationary state for a 
prolonged interval of time. This is true also in the case of 
initial profiles reasonably different from the canonical one. 
In particular, we verified that no (a) = const converges to 
identical slopes as those shown in Fig. A.l at t « 10 12 (cf. 
the bottom row of plots in the Appendix). 

Clear-cut results emerge provided the distribution is 
influenced by a single mode of star-disc interaction. The 
evolution time-scales are shorter in the region where sub- 
stantial fraction of the satellites are already inclined into 
the disc plane, so that rit(a) forgets about the exact ini- 
tial distribution much earlier, at about t <; 10 8 -10 9 . In 
that case the final slope qt ranges from the initial value 
up to the limiting value q^ of the corresponding model and 
the mode of migration (the values are listed in Table [j]) . 
The actual final value of q{ reaches q^ if two conditions 
are satisfied in the region where the power-law distribu- 
tion is established: first, the flow of the satellites in this 
region continues in the disc plane (all bodies are inclined 
to the disc); second, the satellites inflow is faster than the 
changes of their number density at the outer boundary. 
These conditions are satisfied in the late stages until the 
sample is exhausted. In other words, the modified clus- 
ter distribution reflects the adopted disc model while it is 
rather insensitive to the initial form of rio(a). 

This behaviour is illustrated in Figure H in terms of a 
power-law approximation to the computed distributions. 
Given the disc model, different power-law relations may 
be established in separate regions depending on the mode 
of star-disc interaction applicable at that distance. In the 
case (i) and for solar-type satellites, we find that the slope 
qt = —4.9 arises in the region 30^a^l20 (the evolution is 
driven by density waves, accumulating the satellites in an 
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evident peak), while qi — 2.3 suits for 220^aSs680 (i.e. be- 
yond the transition point where the star-disc interaction 
is determined by the gap; here, as an example, qi = 2.3 
is to be compared with the corresponding 04 — y 2 in the 
sixth row in Table 0). Furthermore, the slope is influenced 
by gravitational radiation which dominates the evolution 
on very small radii: qi = 2.8 for r < 20. The three regions 
are distinguished in the graph: gw (dominated by gravi- 
tational waves), dw (density waves), gap (gap formation). 
The large span of values acquired by Of reflects the pres- 
ence of ripples in nt{a) which develop in certain areas of 
the disc. 

The right panel of Fig. @ shows how the rit (a) profile 
is gradually adjusted as time proceeds. Here we plot the 
case (iii) with M* ^ IMq for illustration, and we find 
q = -0.71 (t = 10 10 , 10 < r <, 60), and q = -0.51 
(t = 10 12 , 10 ^ r < 160), respectively. On the vertical 
axis, one can read the fraction of satellites with semi-major 
axes falling in interval (a, a + da). One observes different 
slopes for more massive satellites, M* *J 3Mq, because of 
gap formation (the slope reaches g=0.25); cp. with Fig. || 
to see the effect of M* . On the other hand, the situation 
gets simpler with the other models where such a transition 
does not occur, and a single power-law relation prevails 
at late stages for all M» ^ IOMq. The emerging slopes 
are: (ii) Of = —0.79, and (iv) 0.30, for which we found 
that the quasi-stationary state is established at r ^ 200, 
t = 10 12 . The case (v) develops more slowly and reaches 
the expected value Of = —0.67 somewhat later,[] at time 
t= 10 13 . 



4. Conclusions 

Energy and angular-momentum losses of a satellite col- 
liding with an accretion disc were examined. We explored 
how the relative importance of several secular effects de- 
pends on the parameters of the model, namely, the oscu- 
lating elements of the satellite trajectories, surface den- 
sity of the disc, and the mass and compactness of the 
orbiter bodies. We verified within the thin-disc analysis 
that hydrodynamical drag is, typically, more important 
for the long-term orbital evolution than gravitational radi- 
ation losses. Decay of the orbit due to interaction with the 
gaseous environment brings the orbiter gradually towards 
the center. This drag acts at different regimes depending 
on parameters of the orbiters and of the disc. The magni- 
tude of the drag determines the rate of orbital decay and it 
influences the resulting structure of the cluster, especially 
the characteristic slope of nt{a) density distribution. 

Gravitational radiation losses dominate over the hy- 
drodynamical dissipation if the disc has low density 
and/or the orbiter is of very high compactness (a neutron 
star or a black hole). The two influences are thus comple- 
mentary and they can operate simultaneously at different 



The expected final value of qi for the case (v) coincides 
with qi = —%, as can be seen in the corresponding section 
(gap) of Tab. 



regions of a given source, assuming the thin-disc scheme 
is valid farther away from the center, whereas advection- 
dominated flow takes over at distances of the order of 
w 10 2 r g . One could thus expect the final stages of the 
orbital decay to be governed by emission of gravitational 
waves and by tidal interaction, while the initial transport 
from the outer cluster is ensured by other causes, possi- 
bly the hydrodynamical collisions together with dynamical 
friction. Our results are still only indicative in this part 
because several effects were ignored which must be taken 
into account in realistic models, namely, distribution of 
the satellite orbital parameters may be more complicated 
if two-body collisions between the orbiters are considered. 

Another open question is the impact of star-disc in- 
teractions on the mass function of the satellites. We ob- 
served that M*, E* and £d are the crucial parameters de- 
termining the cluster evolution, but the parameter range 
spans an enormous interval for different types of objects. 
Compact bodies have £* « 10 9 Sq = 1.3 x 10 20 gcm -2 
and hence they are affected less by collisions with the 
disc than solar-type stars, which are aligned with the disc 
plane more rapidly. Consequently, the initial mass func- 
tion is modified towards a higher abundance of compact 
stars residing in inclined orbits, and vice versa. However, 
in the relatively dense environment of the disc, the stars 
accrete at an enhanced rate, they will soon gain sufficient 
mass and eventually collapse, producing additional com- 
pact bodies. A detailed discussion of their subsequent evo- 
lution under the influence of the disc environment remains 
beyond the scope of the present discussion. We only re- 
mark that solar-mass satellites can substantially multiply 
their masses during 10 6 years (Collin & Zahn 1999), as- 
suming that their own radiation and the effect of gaps 
do not halt further accretion. This is almost two orders 
of magnitude shorter than the expected quasar life-time 
(Haehnelt & Rees 1993) and hence the effects of star-disc 
interaction should not be neglected. 

As a concluding remark, let us recall de nouveau that 
the final stages of an orbiter located near the center are 
relevant for gravitational wave experiments. The radia- 
tion losses govern the orbital evolution if the satellite is 
compact enough; on close orbits (e ~ 0.9, a rs 10 2 r g ) the 
influence of gravitational radiation is comparable to the 
effects of star-disc collisions, even if the medium is rela- 
tively dense, Ed ~ 10 5 gcm -2 . 
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Appendix A: Details on modified clusters 

We can expand the results of Sec. |3.2| in slightly more de- 
tail by examining how the modified cluster distribution 
depends on the adopted accretion regime and how the 
quasi-stationary state is eventually established. There is 
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Fig. A.l. A modified distribution nt(a) of satellites in the sample. Columns (i)-(iv) show the modified distributions 
which are developed following the interaction with the corresponding disc models. (The models are listed in Sec. |3|.) 
The last two columns reflect the structure of composite models consisting of different regions, as indicated on top. 
The evolution of the composite models can be compared with the case of individual simple models from which they 
are constructed. In each graph, line types of the curves reflect different views of rit(a) profiles (see the text). First, 
snapshots are shown with varying time (top row; M* = IMq); second, the fraction of satellites aligned with the disc is 
extracted from the whole sample (middle row); third, the quasi-stationary states are plotted for different M„ (bottom 
row). See the Appendix for exact interpretation of different rows of graphs. 



a qualitative difference from the simpler case discussed 
in the main text, because varying surface density profiles 
of the disc models (i)-(v) introduce potentially observ- 
able concentrations of satellites at certain radii. Also, it 
is instructive to recover in the graphs different influences 
resulting from orbital decay, discussed for individual satel- 
lites in previous paragraphs (Sec. H). 



In Figure [A.l] , the first four columns correspond to 
simple models adopted from our list. In the remaining two 
columns, the disc consists of a conjunction of two differ- 
ent models, joined at a given radius: The model denoted 
(ii)+(iii) is formed by the inner (ii) and the outer (iii) part 
of the gas pressure dominated disc with the transition at 
100r g . Similarly, the model denoted (i)+(iv) is a conjunc- 
tion of the radiation pressure dominated disc spreading 
up to 10 3 r g with the outer, self-gravitating disc (Collin 
& Hure 1999). The model parameters determine the final 
form of plotted curves in their mutual interplay together 
with the mass M* of the cluster members, which remains 
the characteristic influencing the efficiency of individual 
contributions to the cluster evolution. 

Given a, the values of nt(a) are distinguished from each 
other if a substantial fraction of the satellites are inclined 
in the plane of the disc without opening a gap in it. In 
such a situation the migration is faster for more massive 



stars whose number density is then lowered with respect 
to low-mass stars. Radial distribution of the satellites is 
influenced by the mode of their migration in the disc, and 
as a consequence the satellites are transported at differ- 
ent rates towards the centre. As a result, a wiggle occurs 
in the rit(a) curve. The particular value of a where the 
satellites are accumulated at an increased rate, forming a 
ring-type structure, depends on the details of the model. 
However, the trend to accumulation at some distance is 
seen in various other situations and it leads to a bias in 
the mass function of the cluster. 



For example, when the model (ii) is considered, the 
stars proceed under the regime of density wave excita- 
tion in the disc plane. The critical point (below which the 
gap would have been opened) lies at too small distances, 
^ 2r g , even for the largest considered mass M* = IOMq. 
Another situation develops in the case (iii) where the gap 
is opened below a — 5, 150 and 7 x 10 3 for our three 
adopted values of satellite masses. This implies different 
shapes of rit{a) for M* = Mq satellites with respect to 
more massive ones. Furthermore, in the marginally unsta- 
ble disc, i.e. the case (iv), the point above which the gap 
could be opened is beyond the outer boundary of the con- 
sidered region and, therefore, the distributions nt{a) are 
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sensitive to the mass of the satellites (most of them are 
inclined in the disc at this stage). 

In order to maintain clarity of the graphs we omit the 
curves corresponding to the model (v). This case is sim- 
ilar to (ii) except for the fact that the gap opens for all 
three considered stellar masses, which makes the case (v) 
distributions very similar to each other, quite independent 
ofM*. 

We could see that simple power-law profiles SdW 
lead to different distributions with respect to more re- 
alistic models, which are governed by different processes 
dominating at the corresponding radius. In discs that are 
a compound of several parts, the surface density profile 
contains several transitions between different regions, and 
hence it is more complicated than a simple power-law. The 
resulting effects on the boundaries resemble the transitions 
which we observed previously with simple discs when the 
mode of radial migration is changed. A local peak of nt(a) 
occurs at radii where the satellites enter the disc region 
with slower radial motion, which is the situation of the 
composite model (ii) + (iii), or quite on the contrary a dip 
develops in nt(a) as it can be seen in the case (i) + (iv). 

We observed that a quasi-stationary distribution n t {a), 
different from the initial one, was established below ;$ 
10 3 r g for rather extended periods of time. However, a de- 
cline of n eventually arrives from large a at final stages 
when the supply of satellites is exhausted. The sample is 



eventually depleted. This can be also seen in Fig. A.l 



the first (top) row, where time evolution is presented for 
one solar-mass satellites. Dimensionless times are t = 10 6 
(solid line), 10 8 (long-dashed), 10 10 (short-dashed), 10 12 
(dotted), and 10 14 (dash-dotted). See footnote 1 for time 
units. The last curve is not visible in all plots because of 
a very large t; it can be noticed e.g. in the case (iv) where 
the evolution is slow and a fraction of satellites persists 
above a<;2x 10 3 . 

Complementary to the above-described panels is the 
second row of graphs where only a part of nt (a) distribu- 
tion is plotted, corresponding to the satellites aligned with 
the disc. Line styles denote time in the same manner as in 
the first row. Naturally, the population of aligned bodies 
is less frequent, and hence some curves are missing in the 
second row with respect to their matching curves in the 
first row, if no satellites are brought to zero inclination. 

Finally, the third (bottom) row shows rit(a) profiles 
taken at t = 10 12 . At this moment a quasi-stationary 
state is already established in the form depending on the 
mode of satellite-disc interaction. The profiles are plotted 
for three different values of the satellite mass, correspond- 
ing to M* = \Mq (thick solid lines), M* = 3M (long- 
dashed lines), and M* = IOMq (short-dashed lines). For 
example, one can check in column (ii) that the structure 
of the plot is ruled by the density wave regime and ac- 
quires the corresponding slope. (The initial slope 1/4 is 
also indicated with a light dotted line.) The curves clearly 
reflect the fact that the orbital decay <M) depends on M* . 
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